function [ pulsesleft, pulsesright, pulsesbase, credperiod ] = periodian( pulsesIn, kThr12 )
%PERIODIAN Summary of this function goes here
%  Detailed explanation goes here

pulsesleft=[];
pulsesright=[];
pulsesbase=[];
credperiod=0;

%% consts
kThr11 = 0.1;
if nargin < 2
    kThr12 = 0.03;
end
kThr2 = 0.03;

%% find out the most probable period
xlast = 0;
pulses = [];
period=0;
noccur=0;
psv_pulses = [];
psv_period=0;
psv_noccur=0;

for i = 1 : length(pulsesIn)
    x = pulsesIn(i);
    if xlast <= 0
        xlast = x;
        continue;
    end
    d = x - xlast;
    if (noccur>0)
        if noccur==1 && abs(d-period)<kThr11*period || abs(d-period)<kThr12*period
            if noccur == 1
                period = d;
                noccur = 2;
            else
                period = period * noccur;
                period = period + d;
                noccur=noccur+1;
                period = period / noccur;
            end
            pulses=[pulses,x];
        else
            if noccur>psv_noccur
                psv_period=period;
                psv_noccur=noccur;
                psv_pulses=pulses;
            end
            period=d;
            noccur=1;
            pulses=[xlast,x];
        end
    else
        period=d;
        noccur=1;
        pulses=[pulses,xlast,x];
    end
    
    xlast = x;
end

if psv_noccur > noccur
    noccur=psv_noccur;
    pulses=psv_pulses;
    period=psv_period;
end

%% get the base
if noccur >= 2
	credbase=pulses(2);
    n = round((pulses(end) - pulses(2)) / period);
	credperiod=round((pulses(end)-pulses(2)) / n);
elseif noccur==1
    credbase=pulse(1);
    credperiod=pulses(2)-pulses(1);
else
    return;
end

%if abs(round(credperiod)-credperiod) > 0.1
%    return;
%end
credperiod = round(credperiod);

pgroupbeginp = pulses(1);
pgroupendp = pulses(end);


%% based on the period just found, accept more pulses that occur at its
%% multiples
for i = 1 : length(pulsesIn)
    x = pulsesIn(i);
    if x < pgroupbeginp 
        m=(credbase-x)/credperiod;
        m = abs(m-round(m));
        if m < kThr2
            pulsesleft=[pulsesleft,x];
        end
    elseif x > pgroupendp
        m=(x-credbase)/credperiod;
        m = abs(m-round(m));
        if m < kThr2
            pulsesright=[pulsesright,x];
        end
    else
        continue;
    end
end
pulsesbase=pulses;
